0.1 Análisis de patrones de puntos

GEMA: CAMBIA LO QUE CONSIDERES, QUE DE TODO ESTO NO TENGO NI IDEA

Las técnicas de análisis de patrones de puntos analizan la distribución de eventos geolocalizados que surgen al azar. La diferencia fundamental con otros análisis que comprenden también el uso de localizaciones (como las temperaturas mínimas medidas por estaciones meteorológicas) es que en este caso los puntos representan eventos conocidos y aleatorios (por ejemplo, los crímenes ocurridos en una ciudad, accidentes de tráfico o incendios en una región). A diferencia de otros eventos, como el ejemplo de las temperaturas mínimas, la ausencia de datos no se debe a la ausencia de medición (e.g. no existe una estación meteorológica en ese lugar), si no a que no se ha producido el evento en dicha localización.

Objetivo de aprendizaje:

El alumno debe ser capaz de conocer los datos de tipo patrones de punto, identificarlos y representarlos adecuadamente.

Tarea 1: Abrimos RStudio

El presente análisis se va a realizar empleando RStudio, por lo que empezaremos abriendo el programa y creando un nuevo script de R en Proyecto>File>New File>R script.

Tarea 2: Importamos y describimos los datos objeto de estudio

El primer paso consiste en importar la base de datos de accidentes de tráfico en la ciudad de Madrid para el año 2020 (Portal de datos abiertos del Ayuntamiento de Madrid). El archivo está en formato csv, por lo que usaremos el paquete readr para importar los datos:

library(readr)
library(dplyr)

# En este caso el archivo está en la carpeta "data" de nuestro proyecto
accidentes2020 <- read_delim("data/2020_Accidentalidad.csv",
  # Configuración adicional
  delim = ";",
  locale = locale(decimal_mark = ",", grouping_mark = "."),
  col_types = cols(
    num_expediente = col_character(),
    fecha = col_date(format = "%d/%m/%Y")
  ),
  na = "NULL"
)


summary(accidentes2020)
#>  num_expediente         fecha                hora          localizacion      
#>  Length:32427       Min.   :2020-01-01   Length:32427      Length:32427      
#>  Class :character   1st Qu.:2020-03-01   Class1:hms        Class :character  
#>  Mode  :character   Median :2020-07-19   Class2:difftime   Mode  :character  
#>                     Mean   :2020-07-07   Mode  :numeric                      
#>                     3rd Qu.:2020-10-22                                       
#>                     Max.   :2020-12-31                                       
#>                                                                              
#>     numero            distrito         tipo_accidente     estado_meteorológico
#>  Length:32427       Length:32427       Length:32427       Length:32427        
#>  Class :character   Class :character   Class :character   Class :character    
#>  Mode  :character   Mode  :character   Mode  :character   Mode  :character    
#>                                                                               
#>                                                                               
#>                                                                               
#>                                                                               
#>  tipo_vehiculo      tipo_persona        rango_edad            sexo          
#>  Length:32427       Length:32427       Length:32427       Length:32427      
#>  Class :character   Class :character   Class :character   Class :character  
#>  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
#>                                                                             
#>                                                                             
#>                                                                             
#>                                                                             
#>    lesividad      coordenada_x_utm coordenada_y_utm  positiva_alcohol  
#>  Min.   : 1.000   Min.   :429069   Min.   :4463495   Length:32427      
#>  1st Qu.: 7.000   1st Qu.:439962   1st Qu.:4471803   Class :character  
#>  Median :14.000   Median :441917   Median :4475026   Mode  :character  
#>  Mean   : 9.862   Mean   :442218   Mean   :4474898                     
#>  3rd Qu.:14.000   3rd Qu.:444328   3rd Qu.:4477652                     
#>  Max.   :77.000   Max.   :454624   Max.   :4488100                     
#>  NA's   :14796                                                         
#>  positiva_droga
#>  Mode:logical  
#>  TRUE:82       
#>  NA's:32345    
#>                
#>                
#>                
#> 

Q1. ¿Qué información contiene nuestros datos?

Este archivo contiene en total 32.427 registros, y proporciona 17 campos asociados a cada registro. Los metadatos de esta información están disponibles en el Portal de datos abiertos del Ayuntamiento de Madrid.

Entre los campos disponibles, podemos destacar los siguientes:

Q2: ¿En qué CRS se encuentran las coordenadas?

Si nos fijamos, las coordenadas no parecen corresponder con longitudes o latitudes, ya que como se explicó el rango posible de valores es \([-180, 180]\) (para longitudes) y \([-90, 90]\) (para latitudes):

Table 1: Accidentes de tráfico en Madrid 2020: Coordenadas
coordenada_x_utm coordenada_y_utm
444597.8 4475156
444597.8 4475156
447710.6 4477156
447710.6 4477156
445076.3 4478372
445076.3 4478372

Tenemos una indicación en el propio nombre de la variable: “UTM” se corresponde con las siglas de “Universal Transverse Mercator,” otra proyección muy utilizada en GIS. En la página web del Ministerio de Agricultura, Pesca y Alimentación podemos comprobar los códigos EPSG habitualmente empleados en la cartografía de España, entre la que se incluyen varias proyecciones UTM.

Como nota, el sistema geodésico de referencia oficial en España es el sistema ETRS891, y el huso UTM de Madrid es el 30 N2, por lo que vamos a probar con el código EPSG 25830, que corresponde con la proyección ETRS89 / UTM zone 30N:


library(sf)
# Objeto sf sin CRS

accidentes2020_sf <- st_as_sf(
  accidentes2020,
  coords = c(
    "coordenada_x_utm",
    "coordenada_y_utm"
  ),
  crs = st_crs(25830)
)

# Comprobamos con la geometría de Madrid

library(mapSpain)
library(ggplot2)
madrid <- esp_get_munic(munic = "^Madrid$") %>%
  st_transform(25830)

# Usamos imagen como mapa de fondo
tile <- esp_getTiles(madrid, "IDErioja", zoommin = 2)

ggplot() +
  layer_spatraster(tile) +
  geom_sf(
    data = accidentes2020_sf,
    col = "blue",
    size = 0.3,
    alpha = 0.3
  )
Accidentes de Tráfico en Madrid (2020)

Figure 1: Accidentes de Tráfico en Madrid (2020)

Q3: ¿Hay algún patrón en la ocurrencia de accidentes de tráfico?

En la figura 1 podemos intuir ciertos patrones en la ocurrencia de accidentes. Por ejemplo, apenas se producen accidentes en la Casa de Campo o en el Monte del Pardo, y parece observarse cierta concentración en las autopistas de salida de la ciudad

Para el siguiente análisis, vamos a analizar el patrón de accidentes en el mes de febrero.

acc_feb <- accidentes2020_sf %>%
  filter(
    fecha >= "2020-02-01",
    fecha < "2020-03-01"
  )

ggplot() +
  layer_spatraster(tile) +
  geom_sf(
    data = accidentes2020_sf,
    col = "blue",
    size = 0.3,
    alpha = 0.3
  )
Accidentes de Tráfico en Madrid (Feb-2020)

Figure 2: Accidentes de Tráfico en Madrid (Feb-2020)

Tarea 3: Análisis de patrones con spatstat

El paquete spatstat ((spatstat_2005?)) es el paquete de referencia cuando se trabaja con patrones de puntos

Además, necesitaremos una ventana de observación (owin). Podemos en este caso obtener los datos espaciales de los distritos de Madrid en la siguiente dirección https://datos.madrid.es/portal/site/egob/menuitem.c05c1f754a33a9fbe4b2e4b284f1a5a0/?vgnextoid=7d6e5eb0d73a7710VgnVCM2000001f4a900aRCRD&vgnextchannel=374512b9ace9f310VgnVCM100000171f5a0aRCRD&vgnextfmt=default

Siguiendo el anterior ejemplo, vamos a analizar el patrón de accidentes en el mes de febrero. Además, en el análisis de patrones de puntos es necesario delimitar la ventana espacial de observación. En este caso será el municipio de Madrid.


library(spatstat)
# Necesitamos un recinto de observación: owin
mad_owin <- as.owin(madrid)

# Seleccionamos del data.frame inicial

acc_feb_df <- accidentes2020 %>%
  filter(
    fecha >= "2020-02-01",
    fecha < "2020-03-01"
  )

# Creamos el objeto ppp
# Es necesario que todas las coordenadas estén proyectadas en el mismo CRS!
feb_ppp <- ppp(
  x = acc_feb_df$coordenada_x_utm,
  y = acc_feb_df$coordenada_y_utm,
  window = mad_owin
)

Q4. ¿Qué información contiene nuestro objeto en formato ppp?

summary(feb_ppp)
#> Planar point pattern:  4047 points
#> Average intensity 6.698035e-06 points per square unit
#> 
#> *Pattern contains duplicated points*
#> 
#> Coordinates are given to 1 decimal place
#> i.e. rounded to the nearest multiple of 0.1 units
#> 
#> Window: polygonal boundary
#> single connected closed polygon with 118 vertices
#> enclosing rectangle: [424891.6, 455899.8] x [4462568, 4499231] units
#>                      (31010 x 36660 units)
#> Window area = 604207000 square units
#> Fraction of frame area: 0.531

Tarea 5: Cálculo de la densidad mediante cuadrantes.

Es importante determinar si los puntos se distribuyen al azar o tienen algún patrón. Por ello, lo primero que haremos será representar el objeto feb_ppp y superponer unos cuadrantes para su comportamiento (véase @ref(fig:plot_ppp_cuadrante)).

## Hallamos los cuadrantes
cuadrante <- quadratcount(feb_ppp,
  nx = 5,
  ny = 5
)
cuadrante
#> tile
#> Tile row 1, col 1 Tile row 1, col 2 Tile row 1, col 3 Tile row 1, col 4 
#>                 0                 0                 0                 0 
#> Tile row 2, col 1 Tile row 2, col 2 Tile row 2, col 3 Tile row 2, col 4 
#>                 0                10                 2                 0 
#> Tile row 2, col 5 Tile row 3, col 1 Tile row 3, col 2 Tile row 3, col 3 
#>                 0                 5                49               730 
#> Tile row 3, col 4 Tile row 3, col 5 Tile row 4, col 1 Tile row 4, col 2 
#>               338                57                 0               163 
#> Tile row 4, col 3 Tile row 4, col 4 Tile row 4, col 5 Tile row 5, col 1 
#>              1488               771                33                 0 
#> Tile row 5, col 2 Tile row 5, col 3 Tile row 5, col 4 Tile row 5, col 5 
#>                73               249                71                 8

## Dibujamos el número de accidentes que hay en cada cuadrante
plot(feb_ppp, pch = 1, main = "")
plot(cuadrante, add = TRUE, col = "red", cex = 0.5)

Como se puede apreciar en la Fig. @ref(fig:plot_ppp_cuadrante)) hay cuadrantes que resgistran cero accidentes y otros qeu registran hasta 1488 accidentes.

Tarea 6: Estimación de la densidad de patrones de puntos.

En la Fig. @ref(fig:intensidad_ripley)) se muestra la gráfica de la estimación usando la K de Ripley. Los conocimientos teóricos necesarios para llevar a cabo este tipo de estimación se verán en el tema Patrones de Puntos. El objetivo de este ejemplo es meramente ilustrativos.

densidad <- density(feb_ppp)
plot(densidad, main = " ")
points(feb_ppp, pch = "+", cex = 0.5)
Intensidad de accidentes de tráfico estimada con la K-Ripley

(#fig:intensidad_ripley)Intensidad de accidentes de tráfico estimada con la K-Ripley


  1. https://www.mitma.gob.es/organos-colegiados/consejo-superior-geografico/csg/etrs89/etrs89-nuevo-sistema-de-referencia-geodesico-oficial-en-espana↩︎

  2. Se puede localizar la zona UTM de una localización mediante la siguiente web: https://mangomap.com/robertyoung/maps/69585/what-utm-zone-am-i-in-#↩︎